%% model_calib.m
%Lee and Maxted; Updated for paper here
%Calibrate model and set up state space
% -------------------------------------------------------------------------

% Numerical assumptions
Nb=20000;
maxit= 100;
crit = 10^(-8);
Delta = 1000;


% Fixed calibration
gamma = 2;
r = 0.03; 
wedgeCC = rCC-r;


% Discount function
beta = 0.75;
rho = 0.05;
betaE = beta;

psiE = (gamma - (1-betaE))/gamma;


% Income process
Nz = 3; %Will not use, since deterministic here

z = 1*ones(1,Nz);
la = 0.001;
la_mat_z = [-la, la/2, la/2; ...
           la/2, -la, la/2; ...
           la/2, la/2, -la];

assert(Nz == 3);
Aswitch = [speye(Nb)*la_mat_z(1,1), speye(Nb)*la_mat_z(1,2), speye(Nb)*la_mat_z(1,3);
           speye(Nb)*la_mat_z(2,1), speye(Nb)*la_mat_z(2,2), speye(Nb)*la_mat_z(2,3);
           speye(Nb)*la_mat_z(3,1), speye(Nb)*la_mat_z(3,2), speye(Nb)*la_mat_z(3,3)];


% State space (construct grid so b=0 is on grid)
bmin = -min(z)/rCC + 1e-12;
bmax = 15;
distanceMinto0 = (-bmin/(bmax-bmin));
b = linspace(bmin,0,round(Nb*distanceMinto0))';
db = b(2)-b(1);
b = [b; db*[1:1:Nb-length(b)]'];
bmax = max(b);

[bb, zz] = ndgrid(b,z);
bb_negative = bb.*(bb<0);

dVf = zeros(Nb,Nz);
dVb = zeros(Nb,Nz);
c = zeros(Nb,Nz);
